

****************************************************************
*** Welfare, alternatives with equivalent energy savings *******
****************************************************************

	**********------------------------------**********
	* Contents: 
		* Compare alternative regulations with same energy savings
		* Benchmarks: flat regulation and actual regulation  
		* Comparison with different values of sigma --> find corresponding kappa 
		* Output: figures displayed in table 6 of paper 

* Load data, based on grid with 4950 cells 

cd "C:\Users\hy65byfe\Desktop\smerge_0712"

* cd "R:\WSV2\TBu_AKe\Spatial_NEW"

use spatial_2_collapse_s, replace    

log using welfare_analysis_table, replace 

keep if year == 2011 // reduce data size 

codebook cell // 4950, but still 50 cells in bunching region.  

sort cell // first bunching (1-50), then restricted (51-1015) then rest. 

gsort -bunching +cell 
gen bcell = _n if bunching == 1 // identify B

gsort -restricted +cell 
gen rcell = _n if restricted == 1 // identify R


*** policy variables ***

scalar E_bar = 1000/220 // benchmark above the max in empirical distribution

scalar K = (51.7*0.59)/220 // scaled to days

scalar kappa = E_bar - `=K' // at compliance line, in days. 

scalar sigma = (47*0.59)/220 // at compliance line, in days


*** parameters: beta, gamma ****
	
scalar beta = 0 // by assumption 

scalar gamma = sigma //  

scalar tolerance = (1/220) // tolerance factor


*** E and EEI in spatial structure ***
gen E_bin = (k_bin/220) + a_bin*`=sigma' // set to midpoint of each cell  
 
replace eei_bin = ((E_bin*220)/(51.7+47*a_bin))* 100 // reset eei_bin to align with compliance line 

*****************************************************
*****************************************************
*** 1) preference parameters


gen a_0 = 		a_bin 						//  pre-ban choice 

gen e_0 = 		E_bar 	- E_bin 			// daily Es, now e from theory. 

gen alpha = 	a_0 	+ beta 	+ e_0*gamma // beta drops out in simplified theory

gen epsilon = 	e_0 	+ a_0*gamma 		//  preference for energy efficiency 


*****************************************************
*****************************************************
*** 2) predicted adjustment 

* preserve 

* gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma)  //  

* gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )	// 

replace restricted = 0 // reset 
replace restricted = 1 if e_0 < kappa - sigma*a_0 - tolerance // restricted = 1 if location is incompliant for given parameters 

gen e_gap = kappa - e_0 - sigma*a_0 // adjustment in e 

gen a_gap = (sigma - gamma)*e_gap    // zero if sigma = gamma 

gen e_star = 	e_gap + e_0 					// post-regulation optimal choice  

gen a_star = 	a_gap + a_0 					// 0 for sigma = gamma.    


sum alpha // preferences in size  

sum epsilon // preferences in energy efficiency  


*****************************************************
*****************************************************
*** 3) evaluate Delta V. 

gen V_gap = 0.5*(e_star - e_0)^2  // formula simplifies when sigma = gamma 

gen E_star = E_bar - e_star  // expressed as consumption


*****************************************************
*****************************************************
*** 4.a) sum up weights (i.e. sales volume)

	*** aggregate 
sort year restricted 

by year, sort: egen w_sum_T =sum(cell_year_sales)  // compute restricted fraction 

by year restricted, sort: egen w_sum_R =sum(cell_year_sales) if restricted== 1 // by year

gen R_ratio = w_sum_R/w_sum_T

by year, sort: egen p_sum_T =sum(cell_year_count) 

by year restricted, sort: egen p_sum_R =sum(cell_year_count) if restricted == 1 // by year
	
gen P_ratio = p_sum_R/p_sum_T
	
*****************************************************
*****************************************************
*** 4.b) sum up Delta V 

	*** at cell level 

gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
sort year restricted 

by year restricted, sort: egen V_sum_R =sum(V_sum_cell) if restricted== 1 // by year

 
sum V_gap if restricted == 1 & year == 2011 // expect: positiv bc V_0 > V_star 

sum V_gap V_sum_cell V_sum_R if restricted == 1 & year == 2011 

*codebook V_sum_R // 14 unique values, one for each year, missings: non-restricted cells.

tab cell a_bin if V_gap < 0 & restricted == 1 & year == 2011 // none.  
 


*****************************************************
*****************************************************
*** 5.a) sum up changes in E  and a 

gen E_gap = E_bin -  E_star // 600 drops out: e_gap = 600 - e_star - 600 + e_0. 

		*** energy savings correspond to "negawatt":  E_gap = E_o - E_star = e_0 + e_star  

sum E_gap a_gap if restricted == 1 & year == 2011 // E_gap must be positive. E_0 > E*. 

gen E_gap_cell = E_gap*cell_year_sales 

sort year restricted 

by year restricted, sort: egen E_sum_R =sum(E_gap_cell) if restricted== 1 // aggregate for year and segment 

gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted consumers 

		*** for a: same process
		
gen a_gap_cell = a_gap*cell_year_sales 
		
sort year restricted 

by year restricted, sort: egen a_sum_R =sum(a_gap_cell) if restricted== 1 // by year, replace 

gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted consumers. 


*****************************************************
*****************************************************
*** 5.b) construct variance 

 	*** at cell level for E 

gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations, multiplied by sales volume in cell. 

	*** aggregate in R 
sort year restricted 

by year restricted, sort: egen D_sum_R =sum(D_sum_cell) if restricted== 1 // sum over R. 

	*** mean of the above. 
gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
sum D_sum_cell D_mean_R if restricted == 1 & year == 2011 // expect: positive when V_0 > V_star 

	*** repeat for a // adjustment is zero for empirical case, < 0 for flat regulations 
	
gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
sort year restricted 

by year restricted, sort: egen DA_sum_R =sum(DA_sum_cell) if restricted== 1 // sum over R. 

	*** mean of the above. 
gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
sum DA_sum_cell DA_mean_R if restricted == 1 & year == 2011 // zero for ABR  


*****************************************************
*****************************************************
*** 5.c) construct coefficient of variation 

 	*** at cell level 

gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


sum cv_R D_mean_R E_mean_R if restricted == 1 & year == 2011 // expect: positive when V_0 > V_star 


*****************************************************
*****************************************************
*** 6) utility ratio  


gen U_ratio = V_sum_R / E_sum_R 

** utility loss / unit of savings ; sum of savings ; sum of loss ; mean loss ; mean variance of loss ; coeffcient of variation; sum of weights ; proportion of restricted product space 

*sum U_ratio E_sum_R V_sum_R E_mean_R D_mean_R cv_R w_sum_R restricted if year == 2011 


sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted if year == 2011 // all outputs reported in table  

local i 1 
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted
	
	
	***** flat ***** 
	
		**** case 1 ****
scalar sigma = 0 // initial: 3.8330
scalar case = 7 // numbering the alternative policy designs 

forvalues i = 3633250(1)3633259 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap 
	
	*replace e_gap = 0 if e_gap < 0 //  
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// 

	gen  a_star = 	a_gap + a_0 					//    

 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  //  
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)

	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  // compute restricted_`i'_`j' fraction 

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate by year and segment 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 // 

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
	*/ 


	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 

	
	gen E_gap = e_star - e_0  

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // E_gap is positive. E_0 > E*. 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 
		
	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations, multiplied by sales volume in cell. 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

	*sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // expect: positive when V_0 > V_star 

	*** repeat for a 
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations
 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 

	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 // output in table 
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}	

* restore  

	*****************************************************
	** 				Simulation over sigma and kappa    **
	*****************************************************
	
		** Note: each case is based on a set sigma, the corresponding kappa is narrowed down iteratively first to arrive at a small range displayed here 
	
		**** case 1 ****
scalar sigma = 0.031 // initial: 3.8330
scalar case = 1

forvalues i = 3832930(1)3832939 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// now what to do with e vs. E 

	gen  a_star = 	a_gap + a_0 					// flawless. tell 'em?   

 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  // 
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)
	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  // 

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 //  

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
 
	*/ 

	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 
	
	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 
		
	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	*** repeat for a // zero for ABR, < 0 for flat regulations 
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 //  
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}

		**** case 2 ****
scalar sigma = 0.063
scalar case = 2

forvalues i = 4033760(1)4033769 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap. 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// 

	gen  a_star = 	a_gap + a_0 					//    

 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  //  
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)

	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  //  

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

 
	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 //  

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
 
	*/ 
 
	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 
	
	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // E_gap must be positive. E_0 > E*. 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 

	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	*** repeat for a // 
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 //  


	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 // 
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}

		**** case 3 ****
scalar sigma = 0.094
scalar case = 3

forvalues i = 4219190(1)4219199 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap. 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// 

	gen  a_star = 	a_gap + a_0 					//    

 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  // careful: general formula needed. 
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)

	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  // 

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

 
	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 //  

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
 
	*/ 

	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 
	
	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // E_gap must be positive. E_0 > E*. 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 

	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	*** repeat for a // zero for ABR, < 0 for flat regulations 
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 // 
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}

		**** case 4 ****
scalar sigma = 0.157
scalar case = 4

forvalues i = 4579390(1)4579399 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap. 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					//

	gen  a_star = 	a_gap + a_0 					//   
 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  //  
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)
	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

 
	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 //  

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
 
	*/ 
 

	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 
	
	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // E_gap must be positive. E_0 > E*. 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 

	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	*** repeat for a //  
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 


	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 //


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 // 
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}


		**** case 5 ****
scalar sigma = 0.188
scalar case = 5

forvalues i = 4742720(1)4742729 { 	
	
	scalar kappa = `i'/1000000
	local k = kappa
	local j = case // first column in table  
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0 - tolerance

	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap. 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// 

	gen  a_star = 	a_gap + a_0 					// 

 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  //  
	
	** test 

	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)
	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  // 

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T

	*by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	

	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

 
	/* 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 //  

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 //  
 
	*/ 


	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 


	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // E_gap must be positive. E_0 > E*. 

	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 

	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year, replace 

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 //  

	*** repeat for a 
	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 

	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 


	*****************************************************
	*****************************************************
	*** 6) utility ratio  
	display kappa 

	gen U_ratio = V_sum_R / E_sum_R 
	
	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 // 
	
	quietly gen w_sum_R_`i' = w_sum_R
	quietly gen E_sum_R_`i' = E_sum_R
	quietly gen E_mean_R_`i' = E_mean_R 
	quietly gen D_mean_R_`i' = D_mean_R
	quietly gen V_sum_R_`i' = V_sum_R
	quietly gen U_ratio_`i' = U_ratio
	quietly gen a_sum_R_`i' = a_sum_R
	quietly gen a_mean_R_`i' = a_mean_R
	quietly gen DA_mean_R_`i' = DA_mean_R
	quietly gen cv_R_`i' = cv_R
	
	quietly gen e_gap_`i' = e_gap 
	quietly gen a_gap_`i' = a_gap
	quietly gen V_gap_`i' = V_gap 
	quietly gen restricted_`i' = restricted_`i'_`j'
* restore 	
	
}

		** weights ** 
gen sw = round(cell_year_sales, 1) // wieghts must be integer
replace sw = 0 if cell_year_sales < 0 // 


save welfare_table_output, replace // intermediate save 

log close 

clear 

use welfare_table_output, replace 

*** checked, get the same values **
log using summary_welfare_table, replace 

*** new values: 
/* 
3633250(1)3633259
3832930(1)3832939
4033760(1)4033769
4219190(1)4219199
4579390(1)4579399
4742720(1)4742729

*/ 
		** flat 	
forvalues i = 3633250(1)3633259 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011 // order of table 01/26/21
} 


		** s = 0.031
forvalues i = 3832930(1)3832939 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011 // order of table 01/26/21
} 

		** s = 0.063
forvalues i = 4033760(1)4033769 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011  // order of table 01/26/21
} 

		** s = 0.094
forvalues i = 4219190(1)4219199 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011  // order of table 01/26/21
} 

		** s = 0.157
forvalues i = 4579390(1)4579399 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011 // order of table 01/26/21
} 

		** s = 0.188
forvalues i = 4742720(1)4742729 {
	scalar kappa = `i'/1000000
	display kappa 
	sum w_sum_R_`i' E_sum_R_`i' E_mean_R_`i' D_mean_R_`i' V_sum_R_`i' U_ratio_`i' a_sum_R_`i' a_mean_R_`i' DA_mean_R_`i' cv_R_`i' if year == 2011  // order of table 01/26/21
} 

*save welfare_table6_2312, replace 

log close

clear 

exit 

